Calculation of the Phase Behavior of Lipids 

M. Miiller* and M. Schick 

Department of Physics, Box 351560, University of Washington 
Seattle, Washington 98195-1560 
*present address: Institut fiir Physik, 
Johannes Gutenberg Universitat, D 55099 Mainz, Germany 

February 1, 2008 



Abstract 

The self-assembly of monoacyl lipids in solution is studied employing 
a model in which the lipid's hydrocarbon tail is described within the Ro- 
tational Isomeric State framework and is attached to a simple hydrophilic 
head. Mean-field theory is employed, and the necessary partition function 
of a single lipid is obtained via a partial enumeration over a large sam- 
ple of molecular conformations. The influence of the lipid architecture 
on the transition between the lamellar and inverted-hexagonal phases is 
calculated, and qualitative agreement with experiment is found. 



1 Introduction 

Lipid bilayers form the framework of biological membranes. Nevertheless al- 
most all membrane lipids adopt non-bilayer configurations, either in their pure 
state or in lipid mixtures, under conditions close to physiological ones ^. 
Although cubic phases are also found ^, the most commonly occurring 
non-bilayer arrangement is the inverted- hexagonalQ, or Hjj, phase in which 
a matrix of hydrocarbon tails is pierced by a hexagonal array of water-filled 
tubes lined by the hydrophilic head groups. Why biological membranes should 
contain lipids which tend to drive it towards a configurational instability, much 
like a lamellar-hexagonal transition, has been the subject of much speculation 
centering on the role such instabilities might play in promoting membrane fu- 
sion, and controlling membrane permeability H, |^. As a consequence of this 
interest, there have been many studies of the phase behavior of lipids in gen- 
eral, and of the parameters which affect the lamellar to hexagonal transition in 
particular||, |^, ||, ^ 0. For example, in a homologous series of saturated dia- 
cyl or diakyl phosphatidylcthanolamincs, an increase of chain length stabilizes 
the Hjj phase, causing the transition temperature between it and the lamellar 
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phase, which exists at lower temperatures, to decreaseQ. Conversely an in- 
crease in the volume of the head group stabilizes the lamellar phase, Lq, causing 
the transition temperature to increase, as clearly demonstrated!^ in mixtures 
of dioleoylphosphatidylethanolamine, DOPE, and dioleoylphosphatidylcholine, 
DOPC. An increase in water content also tends to stabilize Hji. 

A qualitative understanding of these results is provided by a characterization 
of the lipid as a simple geometrical object parameterized by the chain volume, 
the maximum chain length, and the head-group area . The different phases 
result from simple geometrical packing considerations. In contrast to simplifying 
the description of the lipid, large-length scale approaches!^ |l^ simplify the 
description of the bilayer itself, reducing it to an infinitesimally thin membrane 
characterized by elastic constants which are inputs to the theories !Q. Often, 
elements of such theories are combined with other phenomenological terms to 
complete the description !||. 

Simulation of chemically realistic models of lipids and their interactions have 
been carried out|lJ], and yield valuable information about local properties, such 
as density and orientation profiles in the bilayer, lamellar phase, and dynamic 
and transport properties. However simulations are both demanding computa- 
tionally and limited to a rather small number of particles. Within the framework 
of chemically realistic models, the simulation of phase transitions between dif- 
ferent morphologies seems not to be feasible at present. 

Analytic, mean-field, approaches combined with microscopic modeling of 
the tails of lipids, have been applied with success to the manner in which the 
tails pack in the interior of aggregates !^, The contribution of the 

chains to the single-lipid partition function, required by mean-field theory, is 
obtained from an enumeration of the molecular conformations of the tails per- 
mitted within the Rotational Isomeric State (RIS) model!Q. It is here that 
the particular architecture of the chains enters. Bilayer thickness is set by as- 
sumptions on a phenomenological free-energy describing the head-group region. 
These calculations reproduce many features of the density profiles and segment 
orientations in the interior of aggregates. 

The contribution of the lipid tails to the Hn phase has also been examined 
by these means !^1, and it was shown that their entropy always favors this 
phase over the lamellar one. It was also observed that a change in the area 
per head group could lead to a transition to the lamellar phase. No solvent 
was included explicitly, but as its effect would be to alter the area per head 
group, the observation indicates that variation of solvent concentration would 
be able to bring about a transition. The calculation is an approximate one by 
necessity because it is carried out in real space. The full p6mm symmetry of 
the Hji phase is not preserved, as the tubes are considered to be cylinders on 
which the area per head group is uniform. Further, the packing constraints in 
the interstices between the tubes cannot be satisfied exactly. Application of 
this real-space approach to such complicated morphologies as Ia3d would be 
extremely laborious pO|. 
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In this paper, we overcome the difhculties of a real-space approach by com- 
bining the above enumeration techniques with recent advances in the solution 
of mean-field equations of polymers , a combination which we had pre- 

viously shown to be fruitful in a relatively simple systemp3|. Here we apply it 
to a more complicated system, a minimal model of a lipid and solvent mixture, 
one which treats the head group on the same microscopic basis as the tails. We 
are able, thereby, to obtain its phase diagram within mean-field theory and to 
examine how the boundary between and Hu phases changes with lipid ar- 
chitecture, i.e. how that architecture affects the stability of the lamellar phase. 
Our results are in qualitative agreement with experiment. 

2 Model and numerical Self- Consistent Field tech- 
nique 

We describe the lipid as consisting of a tail, comprised of any number of chains, 
(usually two) , containing altogether a total of identical segments of volume 
Vt each, and a head containing two segments of volume Vh/^ each. The solvent 
molecules have volume Vg but are otherwise without structure. The partition 
function of the mixture of Nl lipid molecules and Ns solvent particles in volume 
V can be written 

(1) 

where 

Nl -, 2 

'^" = «^E2E'5(r-r„./.)> (2) 

a=l h=l 
Nl N 

^t = i^s^^<5(r-r„,t), (3) 
are the dimensionless number densities of the head and the tail segments, and 

Ns 

^s^VsY^ S{r - ra,s), (4) 

is the dimensionless number density of solvent particles. The fluid has been 
treated as incompressible. The probability distribution of the lipid configura- 
tions is denoted -P[rQ], and the interaction energy between particles is £. 
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It is convenient to introduce auxiliary fields and consider the particles to in- 
teract with one another via intermediate, fluctuating fields rather than directly; 

J V<^>,-,V<S>tV<i>sVWr,VWtVWsVnexp 1^-^^ (5) 

where the free energy functional, J-, is 

— ■ = A'cln Air In 

keT keT Ns Nl 

-- [ dr \Wh<i>h + Wt<^>t + Ws<i>s +n(l-<i>s~—<^h~-<^>t)}{Q) 

Vs J V Vs J J 

In this expression, Qi^ denotes the single lipid partition function in the ex- 
ternal fields Wh and Wt^ acting on head and tail segments, and Qs the solvent 
partition function; 



[- 2 N 

Ql[Wh, Wt] = I VrP[r] exp <^ -- ^ Wh{rh) - ^ Wt{rt) 

K h—l t—1 



(7) 



Qs [Ws]= J dr exp {-Ws (r) } (8) 
The interaction energy is 

where a and p take the values /i, and s. The incompressibihty condition 

$,(r) + — $^(r) + ^$t(r) = l, (10) 

Vs Vs 

can be used to eliminate the solvent volume fraction in Eq.(^). The terms 
linear in the head and tail volume fractions which result from this procedure 
can be ignored, as they only contribute to the chemical potentials of the heads 
and tails. At this stage we neglect the coupling between the molecular confor- 
mations, the local fluid-like packing of molecular segments, and the local energy 
density and assume that all interactions are contact interactions, 

V^p{r-r')=e^pS{r~r'). (11) 

With these simplifications, the energy can be written in the form 



keT 



^ I |x/it*/i*t - ^Xhh'i'l - ^Xtt^t I , (12) 
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with 
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(13) 
(14) 
(15) 



The functional integration in Eq. (^ cannot be carried out expHcitly. There- 
fore we employ mean-field theory, which approximates the integral by the largest 
value of the integrand. This maximum occurs at values of the fields and densities 
which are determined by extremizing J- with respect to each of its seven argu- 
ments. Fluctuations around these most probable values are neglected. These 
values are denoted below by lower case letters. They satisfy the self-consistent 
equations 



Whir) 
wt{r) 
Ws{r) 
1 

Mr) 
Mr) 
Mr) 



XhtMr) - XhhMr) + ■n{r)vh/vs, 
XhtMr) - XttMr) + TT{v)vtlvs, 
7r(r), 

0s(i') + Mr)vh/vs + Mr)vt/vs, 
NlVs SQl 

Ql Swh{r)' 
NlVs SQl 

Ql Swt{r)' 
Nsvs SQs 

Qs Sws{r)' 



Because the overall density is fixed, we can set j drTr{r) ~ J drws{r) 
mean-field free energy F = J-[4>h, 4>y,4's, Wh, Wt,Ws,TT] is 



(16) 
(17) 
(18) 
(19) 

(20) 
(21) 
(22) 
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XhtMr)(t>tir) + -XhhMr) + 7:Xtt0t (r) 



(23) 



where we have denoted <j)s = NsVg/V, the average, dimensionless, number den- 
sity of solvent, and similarly 4>l = N^Vs/V, for lipids. 

In order to study the self-assembly of lipids into various morphologies, we 
expand the spatial dependence of the densities and fields in a complete set of 
orthonormal functions {/fe(r)}, V^^ J dr fifk — Sik, fi — 1, which possess the 
symmetry of the morphology being considered pl[; e.g. <ph{r) — J2k^h.kfk{r). 
The coefficients 4>h.i, 4>t,ii and i are simply equal to the average, dimen- 
sionless number densities 4>h^ (t>t, £^nd 0^ respectively. The solvent density 
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(?!>s(r) = (j)s,iV expl—Wsir)]/ J dr exp[—Ws]{r) can be Fourier expanded as (/)s,a: = 
^s,i(exp5')fe,i/(exp5)i,i, with 

{S)k,ra = -^W^^j,^ I dr fkfnfm- (24) 
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The partition function of a single Upid in an external field cannot be obtained 
analytically for a realistic architecture. Therefore we approximate the non- 
interacting single lipid probability distribution 7^ by a representative sample of 
Af single lipid conformations. Assigning the Boltzmann weight, tOc, to each lipid 
conformation in the field of mean potentials and wt 



2 Af 



we obtain from Eq. ( ^ ) the components of the tail segment density 

(Pt,k — <Pt,i —jj , (^oj 



E7V 
c=l 



UJr 



and similarly for the head segment density. The self-consistent equations, ex- 
pressed in the basis {fk]^ are solved by a Newton-Raphson-like method. Finally 
we minimize the free energy with respect to the size of the unit cell, D. To this 
end we translate the lipid conformations so as to achieve a uniform distribution 
in the cell. 



3 Results 

The above scheme is applicable to arbitrary lipid architecture and symmetry of 
spatial ordering. We have applied it to model monoolein, a fatty acid whose 
phase behavior in water has attracted much interest |24) . A schematic sketch 
of the architecture is presented in Fig. |^ to illustrate the architectural parameters 
of our model. The single hydrocarbon tail contains N = 17 units, a distance 
i = I.53A apart, with a double bond between the 8th and 9th units. The 
volume of the tail units is Vt — 29 A^. The trans-gauche energy difference 
is taken to be 500cal/mole. We take the two units of the head and the first 
segment of the tail to be collinear, with a distance d between the head units, 
and a distance b between the first tail segment and the adjacent head unit. 
We set all interactions to zero save those between hydrophilic and hydrophobic 
entities, and they are taken to have the same strength e. Thus from Eqs. ( [l^ 
HI), Xhh = 0, xtt = {2evt/kBTvl) = x, and Xht = - Vh)/{2vt). Typically 
between 524,288 and 2,097,152 conformations of single lipids were generated. 
The position of the head was chosen uniformly distributed over the unit cell, with 
a random orientation. The tail was then constructed according to the Rotational 
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Isomeric State (RIS) model||T^. The tail conformations were generated at a 
temperature of 300K, so that they are fluid. 

We use up to 32 basis functions, /fc(r). Because the calculations involve the 
contributions of each individual lipid to the Fourier component of the density, 
the required computer memory scales like the product of lipid conformations 
and number of basis functions. The program needs more than 1 Gbyte memory. 
We employ a massively parallel computer CRAY T3D/T3E and distribute the 
lipid conformations among the processors independently of their spatial position. 
Each processor evaluates the contribution of its assigned lipid conformations to 
the Fourier components of the density according to eqn. (p6|). The partial results 
are collected via shmem routines. The flrst processor also calculates the solvent 
density, and therefore we assign it a smaller number of lipid conformations 
to compensate for the additional work load. Typically between 32 and 128 
processors have been employed in parallel. The program scales linearly with the 
number of processors. 

We first consider the parameters Vt/vg = 1, Vh/vg — 3.2, d/£ — 0, and 
b/£ ~ 1. Thus the head is a relatively small, single, interaction center located 
the same distance from the first tail segment as each tail segment is from its 
neighbors on the chain. The ratio of solvent volume, v^, to lipid volume, vi = 
Vh + Nvt, determines the extent to which the microstructure can be swollen. 
Small solvents tend to swell the microstructure due to their large translational 
entropy, whereas large solvents favor phase separation. The phase diagram 
is determined by calculating the free energy of the different morphologies as 
a function of the solvent concentration and the "temperature" l/x- Phase 
coexistence is determined by equating the chemical potential 



and the Gibbs free energy 



_ dF _ VsdF 
^ ^ dNg ^ Vdct>s' ^ ' 

^-^ (28) 

in the two phases. For a representative value of the temperature, 1/x — 0.25, 
we find the sequence of dimensionless free energy densities, / = VgF/VkgT, 
shown in Fig. ^j. It follows from this sequence that, with increasing solvent 
concentration, there is a transition from the disordered phase to the Hjj phase, 
and from that to the /aSd cubic phase. At larger concentrations, the lamellar 
La phase competes with coexistence between water-rich and laSd phases. This 
is in accord with the phase diagram of monooleinjl^, p^ . 

Having determined that the Ia3d phase does occur in our calculation, we do 
not consider it further. A large number of basis functions is required to deter- 
mine its free energy with sufficient accuracy to determine its phase boundaries. 
Further, it is sufficient to restrict ourselves to the inverted-hexagonal and lamel- 
lar morphologies to investigate the effect of lipid architecture on their relative 
stability, which is our principal interest here. 
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The calculated phase diagram of the system, excluding the cubic phases, is 
shown in Fig. ||. Its salient features are similar to those of the experimental 
monoolein/water mixture [pO|, p^ , and other lipid, solvent mixtures [p5[. The 
Hii phase tends to exist at higher temperatures and water concentrations than 
does the Lq phase. Upon swelling the Hu with solvent at lower temperatures, 
a weak first-order transition to another lyotropic phase, (L^ in this calculation), 
is encountered, while at higher temperatures, a strong first-order transition to a 
disordered, (DIS), solvent-rich phase, is seen. These features are all in agreement 
with experiment. 

The effect of lipid architecture on the L^, Hjj, transition is shown in Fig. ^. 
In (a), we see that at fixed temperature 1/x — 0.25, the solvent concentration 
within the narrow coexistence region between Hn and La phases, shown by 
squares, increases with increasing tail length, while the temperature on the phase 
boundary at fixed concentration, (ps = 0.185, shown by circles, decreases. Thus 
lengthening the tail stabilizes the Hu phase, in agreement with experiment Q. 
The tails in this calculation were taken to be fully saturated so that there 
would be no effect of the relative placement of the double bond. The effect of 
changing the head group volume is shown in (b), and it is seen that increasing 
the volume of the head group destabilizes the Hu phase, again in agreement 
with experiment 

A posteriori, this behavior can also be rationalized in the framework of 
packing models. The morphology is controlled by a packing parameter 77 = 
vi/alc, where a denotes the area per head group and Ic the maximum extension 
of the lipid tail. We set Ic = \/ {R\t) — 15. lA, where {R\i) is the ensemble 
average of the square of the distance between the head unit and the last tail 
segment. The area per lipid tail can be related to the repeat distance D via 
simple geometric considerations: 

2{vh + Nvt) 

aiam = -fTT, : — T lamellar phase (29) 

D{1 - (ps.i) 



Ohex = -fTf. 7 — T \ —7=\(ps + 777 , s mverted-hexagonal phase 

D(l - 0s,i) y V3 V 2{vh + Nvt)J 

(30) 

There is a transition from an inverted-hexagonal to a lamellar phase upon de- 
creasing the packing ratio 77. The molecule is pictured as a wedge with the tails 
constituting the bulkier part. An increase of the head group volume or decrease 
of the hydrocarbon tail length reduces the effective wedge shape of the molecule, 
and therefore tends to stabilize the lamellar phase. 

The occurrence of a transition upon increase of the solvent content is under- 
stood as due to the swelling of the area per head group, which also reduces the 
effective wedge shape of the molecule. The results of our self-consistent field 
calculations are displayed in Fig. ^. Calculating the packing parameter from a 
microscopic model, we find that it decreases upon adding solvent, and takes a 
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value between 0.5 and 0.65 at the transition from the inverted-hexagonal to the 
lamellar phase. From the simple packing arguments, one would have expected 
its value to exceed unity at the transition. This discrepancy only illustrates 
that the phenomenological parameters of the packing model are related but 
approximately to the geometrical parameters of the molecule. 

Along the Hn , transition, we find that the ratio of the lattice constants 
of the two coexisting phases is Dh/D^ ~ I.IO, rather close to the value of 1.16 
extrapolated from the experiments on monoolcin||l^. However the absolute 
values of the lattice spacings are smaller than the experimental ones. The latter 
are = 42A while we obtain 20A. This small value implies that the head 
groups are separated by a very thin water layer, and that hydrocarbon tails 
originating from different monolayers interdigitate significantly, a shortcoming 
encountered in other calculations [p6| . By making the head group bulkier, d/i = 
6, and adjusting the head group volume Vh/vs — 8.5 such that the transition 
occurs at the same solvent density, we increase the calculated result to Dl=28A. 
Extending the head group still further, to take account of their hydration shell, 
would improve the agreement with experiments. 

4 Discussion 

We have explored the self-assembly of monoacyl lipids in solution employing a 
microscopic model. The model describes the architecture of the tail very well, 
by means of the Rotational Isomeric State scheme, and that of the head rather 
crudely. We have calculated the phase diagram employing no other assumption 
than that of mean field theroy. In agreement with experiments, we find that a 
transition from an inverted-hexagonal to a lamellar phase occurs upon increas- 
ing the solvent content or decreasing the temperature, and that an increase 
in the length of the hydrocarbon tail or a decrease in the head group volume 
stabilizes the inverted-hexagonal phase. The ratio of the lattice constants of 
the coexisting phases is in reasonable agreement with experiment. However the 
absolute value of the lamellar spacing is too low. This result highlights one 
of two limitations of our calculation; that we have greatly simplified the head 
group and the solvent. We expect that this deficiency can be remedied with 
more realistic parameterization of the interactions. We have also ignored fluc- 
tuations, which are expected to shift the phase boundaries to somewhat lower 
temperatures and, by inducing effective repulsions between interfaces, to enlarge 
the characteristic length scale of the phases. Because of the extended molecular 
architecture, these effects will be small. They will also have little effect on the 
stability of the lyotropic phases relative to one another. Therefore we are hope- 
ful that our approach, having demonstrated its utility in capturing the effect of 
architecture on lipid phase transitions, will be applicable to the more difficult 
problem of its role in membrane fusion and permeability. 
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Figure 1: Schematic drawing of the Upid architecture of our model. The specific volumes of the head and 
of a single tail segment arc Vh and vt respectively. The distance between the tail segments is b, and that 
between the two segments of the head is d. 
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Figure 2: Free energies at l/xtt = 0.25 of the disordered, La, Ia3d, and Hu phases as function of the 
solvent content. For convenience, the term ^coex4's/kBT has been subtracted from the dimensionless free- 
energy density /, where ficoex is the chemical potential at which solvent-rich and solvent-poor disordered 
phases would coexist were there no stable ordered phases. The dotted line shows the Maxwell construction 
between the and the solvent-rich disordered phase. 
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Figure 4: Influence of lipid architecture on the lamellar to inverted hexagonal transition. The squares 

denote the density 4>s along the coexistence curve at fixed 1/x — 0.25, whereas the circles represent the 
transition temperature, 1/x, at fixed composition (j)s = 0.185. (a) Variation with the chain length N. 
(b) Variation with the head group size Vh- 
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Figure 5: Packing parameter for the lamellar La and inverted- hexagonal Hu phase as a function of the 
solvent concentration for 1/x = 0.25. The vertical line marks solvent composition at which the transition 
occurs. 
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